{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# import necessary modules\n",
    "# uncomment to get plots displayed in notebook\n",
    "%matplotlib inline\n",
    "import matplotlib\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "from classy import Class\n",
    "from scipy.optimize import fsolve\n",
    "from scipy.interpolate import interp1d\n",
    "import math"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# esthetic definitions for the plots\n",
    "font = {'size'   : 16, 'family':'STIXGeneral'}\n",
    "axislabelfontsize='large'\n",
    "matplotlib.rc('font', **font)\n",
    "matplotlib.mathtext.rcParams['legend.fontsize']='medium'\n",
    "plt.rcParams[\"figure.figsize\"] = [8.0,6.0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "common_settings = {'output' : 'tCl',\n",
    "                   # LambdaCDM parameters\n",
    "                   'h':0.6781,\n",
    "                   'omega_b':0.02238280,\n",
    "                   'omega_cdm':0.1201075,\n",
    "                   'A_s':2.100549e-09,\n",
    "                   'n_s':0.9660499,\n",
    "                   'tau_reio':0.05430842,\n",
    "                   'thermodynamics_verbose':1\n",
    "                   }  \n",
    "##############\n",
    "#    \n",
    "# call CLASS\n",
    "#\n",
    "###############\n",
    "M = Class()\n",
    "M.set(common_settings)\n",
    "M.compute()\n",
    "derived = M.get_current_derived_parameters(['tau_rec','conformal_age'])\n",
    "thermo = M.get_thermodynamics()\n",
    "print (thermo.keys())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tau = thermo['conf. time [Mpc]']\n",
    "g = thermo['g [Mpc^-1]']\n",
    "# to make the reionisation peak visible, rescale g by 100 for late times\n",
    "g[:500] *= 100\n",
    "#################\n",
    "#\n",
    "# start plotting\n",
    "#\n",
    "#################\n",
    "#    \n",
    "plt.xlim([1.e2,derived['conformal_age']])\n",
    "plt.xlabel(r'$\\tau \\,\\,\\, \\mathrm{[Mpc]}$')\n",
    "plt.ylabel(r'$\\mathrm{visibility} \\,\\,\\, g \\,\\,\\, [\\mathrm{Mpc}^{-1}]$')\n",
    "plt.axvline(x=derived['tau_rec'],color='k')\n",
    "# The conformal time at reionisation  could be extracted from the code.\n",
    "# But we know it because it is part of the standard output\n",
    "# when thermodynamics_verbose=1\n",
    "plt.axvline(x=4255.316282,color='k')\n",
    "#\n",
    "# Print functions one by one, saving between each (for slides)\n",
    "#\n",
    "plt.semilogx(tau,g,'r',label=r'$\\psi$')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.savefig('thermo.pdf',bbox_inches='tight')"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
